% Empirical exercise 3: portfolio selection


clear all
warning off
clc


% 0. parameter setup
load emp3.mat
TUN = 0:0.01:0.5;                                                                                                                                 
L = [60 120];       % training sizes
c = [1 2];      % contraint value for GEC estimator by Fan et al. (2012)

% 1. forecasting comparison
tic
SAVE = [];
for s = 1:length(DATA)
    DATAt = DATA{s};
    [T,N] = size(DATAt);
    DATAt = DATAt - repmat(RF,1,N);   % substract risk-free rate
    for i = 1:length(L)
        WL = L(i);
        SR = [];
        parfor t = 1:(T-WL)/12   % roll one-year at a time
            int = 12*(t-1)+1:12*(t-1)+WL;
            ine = 12*(t-1)+WL+1:12*t+WL;
            Xt0 = DATAt(int,:);
            Xe = DATAt(ine,:);         
            Xt = Xt0-repmat(mean(Xt0,1),size(Xt0,1),1);       % de-mean here
            % 1.1 l2-relax 
            W = zeros(N,3);
            for opt = 1:3
                W(:,opt) = l2relax_ps(Xt, TUN, opt);
            end
            % 1.2 GEC by Fan et al. (2012)
            W2 = zeros(N,length(c));
            for j = 1:length(c)
                W2(:,j) = gec_fun(Xt, c(j));
            end
            % 1.3 l2-relax without summing to 1 restriction
            W3 = zeros(N,3);
            for opt = 1:3
                W3(:,opt) = l2relax_ps2(Xt, TUN, opt);
            end
            
            % 2. construct various portfolios
            P0 = mean(Xe,2);            % simple averaging
            P1 = Xe*W;                  % l2relax (opt = 1,2,3)
            P2 = Xe*W2;                 % GEC
            P3 = Xe*W2;                 % l2relax no summing to 1 restriction

            % 3. evaluate by Sharpe Ratio
            P = [P0 P2 P1 P3];
            SR(t,:) = mean(P)./std(P);
        end
        SAVE = [SAVE; mean(SR)];
    end    
end
toc

% 4. show results
VN = {'60','120'};
prt_tab(SAVE,repmat(VN,1,3),3)

% 5. plot weights
for j = 1:3
    DATAt = DATA{j};
    [T,N] = size(DATAt);
    DATAt = DATAt - repmat(RF,1,N);   % substract risk-free rate
    L = [60 120];       % training sizes
    TUN = 1;
    IN = [0.15 0.05 0 -0.05];
    VN = {'(a) Five-year-training','(b) Ten-year-training'};
    figure (2+j)
    for i = 1:length(L)
        Xt = DATAt(1:L(i),:);    
        w = l2relax_ps(Xt, TUN, 1);
        subplot(1,2,i)
        plotline(IN,sort(w,'descend'),'Portfolio')
        title(VN{i})
        if j==1
            ylim([-0.2501 0.2501])
        else
            ylim([-0.1502 0.1502])
        end
    end
end

%% nested function
function w = l2relax_ps2(X, TUN, opt)
% L2-relaxation estimation for portfolio selection without summing to 1 restriction. 
% (nest CV by the Sharpe ratio)
%
% TUN = tunning parameter (can be a scalar or vector)
% opt = 1. no shrinkage
%       2. linear shrinkage
%       3. nonlinear shrinkage

    % 1. nest CV or one-time estimation
    if length(TUN) > 1  
        SR = zeros(length(TUN),1);
        for i = 1:length(TUN)            
            SR(i,1) = nest_cv_ps(X,TUN(i),opt); % nest CV with Sharpe ratio
        end        
        % 2. find the best
        [~,IN] = max(SR);       % choose the largest sharpe ratio
        tun = TUN(IN);
    else
        tun = TUN;
    end
    
    % 2. estimate with optimal tunning
    w = l2relax2([],X,tun,opt);        
end

function RISK = nest_cv_ps(X,tun,opt)
    % 1. construct 5 nested index
    n = size(X,1);
    d = round(n/5);
    tmp = 1:d:n;
    IN = [tmp(1:5)' [tmp(2:5)'-1; n]];
    INt = [ones(4,1) IN(1:4,2)];
    INe = IN(2:end,:);
    if d <= 12         % to avoid the "analytical_shrinkage.m" error
        INt = INt(2:end,:);
        INe = INe(2:end,:);
    end
    % 2. obtain risk for each tuning value
    RISKt = zeros(size(INt,1),1);
    for j = 1:size(INt,1)
        Xt = X(INt(j,1):INt(j,2),:);        
        Xe = X(INe(j,1):INe(j,2),:);
        wt = l2relax2([],Xt,tun,opt);     
        % 2.1 estimate (rough) Sharpe Ratio 
        r = Xe*wt;
        RISKt(j) = mean(r)/std(r); 
    end
    RISK = mean(RISKt);    
end
 

